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Abstract  Given  a triangulation  of  a planar  region,  the  reduced  Hsieh-Clough-Tocher  (rHCT) 
triangular  element  enables  the  construction  of  a smooth  piecewise-cubic  surface.  In  prepa- 
ration for  the  use  of  energy  minimization  to  select  interpolants  from  that  family  of  surfaces, 
a formula  for  the  bending  energy  of  the  generic  rHCT  element  is  determined,  treating  the 
element  as  a thin  and  almost  flat  plate.  The  aim  is  to  find  an  approach  to  interpolation  which 
displays  low  sensitivity  to  the  choice  of  triangulation.  The  problem  arose  in  terrain  modeling. 

Keywords:  Delaunay  triangulation,  finite  elements,  Hsieh-Clough-Tocher  element,  interpo- 
lation, spline,  surface,  terrain  modeling,  triangulation 


1.  Introduction 

Lawson  [9], [10]  pioneered  the  use  of  (^-compatible  finite  elements  for  the  task  of  interpolating 
a continuously  differentiable  (“C1”)  function  z — f{x,y)  from  a finite  set  of  spatial  points 

(x i,Vi,Zi),  i = 1 

so  that  Zi  = z(xi,yi).  An  important  application  is  to  terrain  modeling  from  large  sets  of 
elevation  data,  where  the  function  values  Zi  are  elevations  at  specified  data  locations  (xl,yl). 

For  his  work,  Lawson  chose  an  element  generally  ascribed  to  Clough  and  Tocher  [4],  but 
frequently  referred  to  in  the  literature  (e.g.  Ciarlet  [2])  as  the  reduced  Hsieh-Clough-Tocher 
(rHCT)  element  For  a description  of  that  and  related  elements,  we  recommend  the  text 
by  Bernadou  and  Boisserie  [1].  Farin  [7]  formulated  Clough-Tocher  interpolation  in  terms  of 
Bernstein-Bezier  polynomials  and  also  provided  a modification  of  the  rHTC  element  which 
trades  C2-continuity  at  the  barycenter  - where  subtriangles  with  different  cubics  join  (see 
Sections  2 and  4 below)  - for  improved  approximation  to  C2-continuity  at  the  edges  of  the 
triangular  element.  For  general  information,  the  reader  may  want  to  consult  Ciarlet  [3]  and 
Zienkiewicz  [18],  to  name  just  two  representatives  of  a large  body  of  literature  on  the  Finite 
Element  Method. 

The  rHCT  element  is  a triangular  surface  patch,  that  is,  it  is  defined  over  a footprint 
triangle  in  the  x,  y-plane.  Its  use  'in  the  context  of  interpolation  thus  presupposes  a “ trian - 
gulation ” of  the  data  locations  (x^,?/*),  i = l,...,n,  that  is,  a set  of  triangles  whose  interiors 
do  not  meet,  whose  vertices  are  the  data  locations,  and  whose  union  covers  the  convex  hull 
of  the  latter.  Two  different  triangles  of  the  triangulation  thus  have  either  a single  edge  or  a 
single  vertex  in  common,  or  do  not  meet  at  all.  In  the  context  of  terrain  modeling,  a set  of 
elevation  data  locations  triangulated  in  this  fashion  is  now  frequently  called  a “ triangulated 
irregular  network  (TIN)”.  Any  set  of  planar  data  locations  (. Xi,yi ) can  be  triangulated  in 
many  different  ways  and  this  affects  the  surface  to  be  constructed.  A frequently  used  generic 
method  is  the  Delaunay  triangulation  [5]. 
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The  rHCT  elements  - by  their  construction  - join  smoothly  at  triangle  boundaries,  thus 
ensuring  a continously  differentiable  fit.  The  specification  of  an  rHCT  element  over  a particular 
triangle  requires  that  the  function  value  Zi  and  the  gradient  (zjX,  Ziy)  be  given  at  each  triangle 
vertex,  that  is,  at  each  data  location  ( x^yi ).  The  resulting  function  is  to  agree  with  those 
prescribed  values  and  slopes. 

Unless  such  gradient  information  is  also  provided,  it  must  be  estimated  from  given  function 
values  Zi  in  order  to  complete  the  specification  of  rHCT  elements  for  each  triangle  of  the  TIN. 
Thus  one  distinction  between  alternate  methods  for  surface  interpolation  based  on  the  rHCT 
element  is  what  method  for  estimating  gradients  at  data  locations  is  chosen. 

Lawson’s  approach  [10]  is  to  select  at  least  six  neighbors  of  a data  point,  assign  them 
weights  which  decrease  with  distance  from  the  data  point,  and  then  consider  6-parameter 
quadratic  bivariate  functions  which  pass  through  the  data  point.  From  among  those  functions 
select  the  one  which  minimizes  the  weighted  least  squares  error.  The  tangent  to  that  quadratic 
function  at  the  data  point  provides  the  gradient  estimate.  Franke  [8]  and  Stead  [15]  compare 
various  procedures  for  triangulation  based  interpolation,  in  particular,  aspects  of  gradient 
estimation. 

In  the  context  of  terrain  modeling,  with  elevation  data  given  along  digitized  contour  lines, 
Mandel,  Bernal,  and  Witzgall  [11]  arrived  at  gradient  estimates  by  minimizing  the  elastic 
energy  of  a mechanical  surrogate  structure  for  the  surface  consisting  of  thin  beams  along  the 
edges  of  the  triangular  patches  and  joined  tangentially  to  thin  plates,  one  at  each  vertex. 
The  orientation  of  the  plates  can  be  adjusted  to  minimize  the  elastic  energy  of  the  beams, 
thus  defining  gradients  at  triangulation  vertices,  that  is,  at  data  locations.  That  procedure 
was  chosen  because  of  the  distribution  of  data  points  along  lines,  which  rendered  local  fitting 
procedures  less  attractive.  It  also  eliminated  the  somewhat  arbitrary  choice  of  neighboring 
data  points  and  their  weights. 

The  two  rHCT  methods  mentioned  above  are  less  sensitive  to  changes  in  the  underlying 
triangulation  than  the  still  most  frequently  employed  linear  TIN  methods.  In  those  linear 
methods,  each  footprint  triangle  in  the  TIN  gives  rise  to  a planar  triangular  facet  spanned 
by  the  elevations  at  the  triangle  vertices.  The  resulting  piecewise  linear  surface  is,  indeed, 
extremely  sensitive  to  the  choice  of  triangulation. 

With  the  goal  of  further  reducing  the  sensitivity  of  the  interpolating  surface  to  the  choice 
of  triangulation,  we  propose  to  examine  an  alternate  method  for  estimating  the  gradients  at 
data  locations.  The  idea  is  to  replace  the  surrogate  mechanical  structure  consisting  of  thin 
beams  used  in  [11]  by  a surrogate  mechanical  structure  consisting  of  the  actual  rHCT  elements, 
joined  together  smoothly,  and  to  minimize  the  total  elastic  energy  of  the  resulting  interpolating 
surface  considered  as  consisting  of  almost  flat  thin  plates.  More  precisely,  we  determine  that 
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interpolating  surface  z — z(x,y)  which  consists  of  rHCT  elements  and  minimizes 


(i.i) 


d2z 
dx 2 


+ 2 


a2. 


dx  dy 


+ 


d2. 


dy'4 


dx  dy 


with  respect  to  the  choice  of  gradients  at  triangulation  vertices. 

A critical  first  step  in  this  direction  is  to  develop  closed  formulas  for  that  energy  operator 
as  applied  to  the  generic  rHCT  element,  and  that  is  the  purpose  of  this  report.  Those  formulas 
are  unwieldy,  and  the  availability  of  symbolic  computation  packages  was  instrumental  in  their 
derivation.  This  work  relied  on  Mathematica  (see  for  instance  [16]). 

Powell-Sabin  splines  [14]  also  provide  a unique  piecewise  polynomial  Ci-function  which 
meets  prescribed  values  and  gradients  at  the  vertices  of  a triangulation.  The  concept  of  thin 
plate  energy  minimization  has  also  been  considered  for  various  other  approaches  to  bivariate 
interpolation  (e.g.  Powell  [13]).  Mansfield  [12]  minimizes  the  full  strain  energy  functional. 
The  reader  may  want  to  consult  Dierckx  [6]. 

We  also  revisit  the  definition  of  the  rHCT  element  and  verify  its  main  properties.  This 
exposition  as  well  as  the  derivation  of  the  energy  formulas  is  conducted  in  terms  of  barycentric 
coordinates.  In  this  fashion,  the  inherent  symmetries  of  the  element  are  preserved. 

In  particular,  the  roles  of  the  vertices  of  the  element  and  their  associated  quantities  are 
interchangeable.  More  precisely,  many  formulas  involving  indexed  quantities  follow  from  each 
other  by 

(1.2)  cyclic  substitution  : 1 — 2 — ^ 3 — > 1, 

where  each  index  value  is  replaced  by  its  cyclic  successor. 

We  will  not  always  write  all  three  instances  of  formulas  that  arise  from  each  other  by  cyclic 
substitution  of  indices.  Instead,  we  will  write  one  instance  of  the  formula,  and  indicate  that 
the  remaining  instances  can  be  generated  by  cyclic  substitution. 


2.  Definition  of  the  rHCT  Element 

In  this  section  we  define  various  quantities  which  are  used  to  calculate  the  generic  element. 
Suppose  we  are  given  three  triangle  vertices, 

(£i,2/i),  (£2,2/2),  (£3,2/3), 

with  associated  function  values,  i.e.,  elevations, 


*i,  z2,  z3, 

and  partial  derivatives  with  respect  to  x and  y,  i.e.  coordinate  slopes, 

Zlx,  ^lyi  ^2xi  %2i /,  xi  Z3y. 
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We  wish  to  define  the  elevation  z — z(x,  y ) at  any  location  (x,  y)  in  the  footprint  triangle  of 
the  rHCT  element,  so  that 

dz  dz 

Zj  Z,{Xi,  yi),  Zix  (^Xx  5 yi)  i Ziy  (*^z : Vi)  ■>  ^ 1,2,3. 

For  that  purpose,  it  is  convenient  to  express  functions  over  triangles  using  ubary centric 
coordinates” . For  any  point  (x,  y)  in  the  plane,  its  barycentric  coordinates 


Ai,  A2,  A3 


are  defined  by  the  relations 

Ai  + A2  4-  A3  = 1 

(2.1)  X\X\  + A2X2  + A3X3  = x 

A12/1  + A22/2  + A3y3  = y . 

Following  Zienkiewicz  [17],  we  use  the  abbreviations 

Xij  :=  Xi  - xj:  yio  :=  yx  - yjx  % :=  z{  - Zj , ij  = 1,  2,3,  i / j. 

The  following  determinants  will  play  a role: 


(2.2) 


'xy 


zy 


1 

1 

1 

a 

Q 

I 

II 

:=  det 

Xi 

x2 

23 

V\ 

2/2 

2/3 

1 

1 

1 

Q 

I 

II 

:=  det 

Z\ 

22 

23 

V\ 

2/2 

2/3 

= -Dxz 

:=  det 

1 

1 

1 

2l 

22 

23 

Xi 

x2 

£3 

— 242/23  + ^22/31  + x3Vl2, 


Z\y23  + ^2  2/31  + 232/12: 


Z\X^3  + 22^31  + 23X12. 


The  determinant  Dxy  is  double  the  area  of  the  footprint  triangle.  It  is  assumed  that  the  area 
of  that  triangle  does  not  vanish.  The  footprint  triangle  is  required,  moreover,  to  have  positive 
orientation,  that  is, 


Dxy  > 0. 
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Solving  the  linear  system  of  equations  (2.1)  yields 


(2.3)  Xi  = {y23X  + x32y  + x2ys-y2X3)/DXy 

X2  — {y3ix  + x\sy  + x^yi  — yzX\)IDxy 
A3  — (ynx  + x2iy  + x\y2  ~ yix2)/Dxy. 


Note  that  the  above  formulas  arise  from  each  other  by  cyclic  substitution  (1.2).  Indeed, 
barycentric  coordinates  treat  all  vertices  of  a triangle  symmetrically.  Their  signs  indicate 
whether  a location  lies  inside  or  outside  the  footprint  triangle:  a negative  barycentric  coor- 
dinate means  the  location  is  outside  the  triangle;  if  all  three  coordinates  are  positive,  the 
location  lies  in  the  interior  of  the  triangle. 

Edges  of  the  footprint  triangle  are  characterized  by  the  vanishing  of  one  of  the  barycentric 
coordinates.  At  vertex  i,  A*  = 1 and  X0  = 0 for  j / i.  The  barycenter  or  centroid 


(xo,  2/o ) 


f x\  + x2  + ^3  Vi  + 2/2  + yz 
V 3 5 3 


of  the  triangle  is  characterized  by  Ax  = A2  = A3  = 1/3. 

The  barycenter  plays  a key  role  in  the  definition  of  the  rHCT  element.  It  is  used  to  define 
a “ barycentric  partition ” of  the  triangle  into  four  subsets,  the  barycenter  itself  and  three 
subtriangles  (see  Figure  1): 


Bo  '=  { (Ai , A2,  A3)  : Ai  = A2  = A3  = 1/3} 

B\  :=  {(Ai,A2,A3)  : 0 < Ax  < A2,  Ai  < A3} 

B2  :=  { (Ai,  A2,  A3)  : 0 < A2  < A3,  A2  < Ai } 

B3  { (Ai,  A2,  A3)  : 0 < A3  < Al5  A3  < A2}  . 


Note  that  in  each  of  the  subtriangles  Bi,i  = 1,  2,  3,  the  corresponding  barycentric  coordinate 
is  dominated  by  the  remaining  ones;  e.g., 

(2.4)  Ax  < min{A2,  A3}  in  B\. 

The  rHCT  function  is  defined  in  terms  of  three  correction  functions  [10] 


Pi,  P2,  P3, 

which  are  piecewise  cubic  polynomials  in  the  barycentric  coordinates,  each  defined  in  a piece- 
wise  manner  with  respect  to  the  subtriangles  of  the  barycentric  partition. 
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81 


(2.5) 


Pi  •= 


for  (Ai,  A2,  A3)  € Bq 


A1A2A3  + -A*  — — Aj  for  (Ai,  A2,  A3)  € B\ 
0 2 


—A  l + 1a?a 


2^3 


— -A?  + -AoA 
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3^2 


for  (Ai,  A2,  A3)  € B2 


for  (Ai,  A2,  A3)  € B$. 


Expressions  for  the  remaining  correction  functions  p2  and  p3  follow  directly  from  the  above 
expression  by  cyclic  substitution  (1.2). 

The  edge  of  the  footprint  triangle  spanned  by  vertices  1 and  2 is  characterized  by  A3  — 0, 
and  essentially  lies  in  subtriangle  B3.  Since  the  Z?3-expression  for  each  correction  function 
contains  A3  as  a common  factor,  those  functions  vanish  along  that  edge.  In  general, 


(2.6) 


Pi  = P2  = p3  = 0 on  the  triangle  boundary. 


Each  function  p*,  z = 1,2,3,  is  continuous  across  the  subtriangle  boundaries.  This  is  readily 
verified  for  the  correction  function  px  as  follows.  At  the  boundary  between  subtriangles  J52  and 
B$,  we  have  A2  = A3.  If  we  substitute  the  single  quantity  A for  these  values,  the  expressions 
for  pi  in  B2  and  B$  become  identical.  At  the  boundary  between  subtriangles  Bi  and  B2,  we 
have  similarly  Ai  = A2  :=  A.  Substituting  A into  the  expressions  keyed  to  B\  and  B2  yields 

-\2  - 1\3 
2 6 

in  both  cases.  The  argument  is  analogous  for  the  remaining  two  correction  functions.  It 
will  be  seen  in  the  next  section  that  the  correction  functions  are  also  smooth  at  subtriangle 
boundaries,  and  thus  belong  to  the  class  C1. 

The  following  quantities  can  be  calculated  easily  from  the  initial  triangle  data. 


(2.7) 


Mji  ZixXji  T Ziyl/ji,  i,  j — 1,  2, 3,  i ^ j ■) 


Mji  + Mij 


Ci: j = 


Qij  2 7 ~ 2 

Li  = \Jx)k  + y% > 2 # j , j # k + i, 

3 m-Ll)  „ 3 (L\-L\) 


— z 


jz 


K i = 


L\ 


Ko  = 


L\ 


K,  = 


3(1?  - L\) 


L\ 


The  six  quantities  Mji  represent  the  directional  derivatives  at  vertex  i with  respect  to 
the  direction  vector  (x^-,  yZJ),  which  represents  an  entire  directed  edge  of  the  triangle.  Note 


that  for  a linear  function  z — z(x,y),  Mji  — —Mio  = Zij , and  for  a quadratic  function, 
Mji  — Mij  = 2 Zji.  As  a consequence,  the  quantities  Qi3  and  Ci3  both  vanish  for  linear 
functions,  and  the  quantities  Ci3  for  quadratic  functions.  Also  note 


Qji  — Qij->  and  Cji 


We  also  define  three  functions 


Vi,  v2,  v3 

in  terms  of  the  previously  (2.5)  introduced  correction  functions  p1?  p2>  P3: 

(2.8)  V\  :=  A2A3(A2  — A3)  + K\P\  — p2  + p3 

Expressions  for  V2,  V3  follow  by  cyclic  substitution  (1.2).  The  functions  Vi  are  continuous 
and,  by  (2.6), 

(2.9)  V\  — V2  = V3  = 0 at  each  triangle  vertex. 

With  the  preceding  definitions,  we  are  now  ready  to  present  a formula  for  the  rHCT 
function,  expressing  the  elevation  z = z(x,  y ) at  any  location  (x,  y ) in  the  footprint  triangle: 

z — Z\\\  22A2  + z3 A3  + 

(2.10)  Q23  A2  A3  + Q3lA3Ai  + Ql2AiA2  + 

C23V1  + C31V2  + C12 V3  . 

This  indeed  defines  a function  z = z(x,y ),  since  the  quantities  A * are  functions  of  x and  y 
by  (2.3).  To  see  that  z = z(x,y)  assumes  the  prescribed  elevations  Z{  at  the  vertices  of  the 
footprint  triangle,  recall  (2.9)  and -that,  at  vertices,  one  A;  equals  1 while  the  others  are  zero. 

Recall  also  that  for  a quadratic  function  the  coefficients  CZJ  vanish,  so  that  the  first  two  lines 
of  the  above  expression  for  z — z(x,y)  describe  an  interpolation  that  is  quadratic.  Similarly, 
if  both  sets  of  coefficients  Cij  vanish,  then  the  first  three  terms  of  (2.10)  describe  a 
linear  interpolation.  This  observation  is  significant  because  the  energy  integral  (1.1)  vanishes 
for  linear  functions,  and  will  thus  be  a homogeneous  quadratic  form  in  the  six  coefficients 

Qij ■>  Cij. 

3.  First  Derivatives  of  the  rHCT  element 

In  this  section,  the  partial  derivatives  of  2 = z(: r,  y)  will  be  determined,  enabling  us  to  establish 
continuous  differentiablility.  The  derivatives  of  the  rHCT  function  z — z(x,  y)  are  in  the  final 
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instance  based  on  the  derivatives  of  the  barycentric  coordinates  A *. 


(3.1) 


d\i 

dx 

dX2 

dx 

8X3 

dx 


2/23 

dXi 

^23 

D 

^xy 

dy 

D 5 

2/31 

dX2 

X31 

B xy 

dy 

D ’ 

yx 

2/21 

d\3 

X21 

D ’ 

^ xy 

dy 

D yx 

:y  and  D 

yx 

-Dxy  as  defined 

permits  us  to  write  the  above  formulas  symmetrically  in  x and  y. 

For  the  correction  functions  pi  defined  in  (2.5),  the  Chain  Rule  yields 


and 


(3.2) 


(3.3) 


dpi 

dx 


Similarly, 


dpi 

dy 


dpi 

dx 


18 


dpi  dXi  dpi  dX2  dpi  d A3 
dXi  dx  + dX2  dx  + dX3  dx 


y 23 


D 


xy 


^2^3  + ~X\  — Ai  ) 2/23  + A1A3  2/31  + A1A22/12 
+ ^2 A3  ) 2/31  + ~ X2  yi2 


7^3  7/31  + -A3  + A3A2^  2/12 


1 1 
^23j 


18 


D 


yx 


^A2A3  4-  ~-X\  — Ai^  x2s  4-  A1A3  X31  4-  AiA2Xi2 
(^—  7^2  + A2A3^  X31  4-  ~X\xi2 
-A3 x31  + -Ag  + A3A2j  x12 


in  B0 
in  Bi 
in  B2 
in  B3 


in  B0 
in  Bi 
in  B2 
in  £3 


The  above  two  formulas  are  symmetric  in  x and  y.  The  corresponding  formulas  for  the 
derivatives  of  the  other  correction  functions  p2  and  p3  follow  by  cyclic  substitution  (1.2). 
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It  is  also  readily  seen  that  the  gradient  of  pi  vanishes  along  the  R2-edge  - characterized 
by  A2  = 0 - and  the  B3-e dge  - characterized  by  A3  = 0 - of  the  footprint  triangle  (compare 
2.6): 

(3.4)  ^ = d-p-  = 0 for  A;  = 0,  i = 1,2,3. 

dx  dy 

We  are  now  able  to  ascertain  that  the  correction  functions  pi  are  indeed  continuously 
differentiable,  that  is,  that  they  join  smoothly  at  subtriangle  boundaries.  We  consider  only 
pi.  The  other  cases  follow  by  cyclic  substitution  (1.2).  The  boundary  between  subtriangles 
B2  and  Bs  is  characterized  by  A2  = A3.  The  expressions  specified  in  both  (3.2)  and  (3.3)  for 
those  subtriangles  are  identical  in  this  case.  At  the  boundary  between  Bi  and  B2  we  have 
Ai  = A2  :=  A and  A3  = 1 — 2A.  Applying  these  substitutions  for  subtriangle  B i yields 

(A(l  — 2A)  + b-A2  — A)y23  4-  A(1  — 2A )?/31  + X2y\2 


or,  taking  into  account  that  y23  = — yi2  — y31, 


(— 2^2  + ^)2/3i  + 2 ^2yi2 ■ 


The  same  expression  clearly  results  in  subtriangle  B2.  The  argument  for  the  derivative  with 
respect  to  y is  analogous,  and  so  is  the  consideration  of  the  boundary  between  B\  and  £3. 

Since  the  correction  functions  pi  are  continuously  differentable,  the  same  holds  for  the 
functions  V*  in  formula  (2.10)  for  the  function  2 = z(x,y).  As  a consequence,  the  rHCT 
element  belongs  to  the  class  Cl. 

Using  the  Chain  Rule  and  formulas  (3.1),  we  find 


dz_ 

dx 


Thus, 

(3.5) 


dz 

dx 


,i  dz  dX 
- + 


,2  dz  dX c 
dX2  dx  + 


<9A3  dx 
l + Q3i^3  + QnX2)y23  + 


dV> 


+ (z3  + Q23X2  + Q3lX\)yi2]  + 


~jZ. — [^1^/23  + z2y3l  + z3y\2  + 

Xi{Qny3i  + Q312/12)  + X2(Q23yu  + Q\2y2s)  + A3(Q3i?/23  + Q23?/3i)]  + 


C 


23 


+ C31 


dx 


+ C\2 


dx 
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Exchanging  x and  y then  yields 

(3.6)  — = — — [Z\X2Z  + ^2^31  + 2:3X12  4- 

dy  DyX 

^l(0l2^13  + 031^2l)  + ^2(023^21  + Q 12*^32)  + ^3(^31^31  + Q232/13)]  + 


G 


23 


dV 1 
dy 


+ C ; 


31 


dV2 

dy 


+ C*i2 


m 

dy 


Note  that  the  constants  in  both  expressions  equal  the  determinants  Dzy  and  Dzx , respectively, 
as  defined  in  (2.2). 

The  derivatives  of  the  functions  Vi,  i — 1,2,3,  (2.8)  are  next.  Again,  we  display  only  the 
derivatives  of  Vi,  since  the  expressions  for  the  derivatives  of  the  remaining  functions  follow  by 
cyclic  substitution  (1.2).  The  first  term  of  the  function  V\  is  the  polynomial  A2A3(A2  — A3). 
By  (3.1),  the  partial  derivative  of  this  expression  with  respect  to  x is  given  by 

— — [(2A2A3  — A fjysi  + ( A|  — 2A2A3)yi2]. 

-L^xy 


Thus,  - and  symmetrically  for  the  derivative  with  respect  to  y -, 


(3.7) 


dV 1 
dx 
dV1 
dy 


— — [(2A2A3  — A3)?/3i  + (Ag  — 2A2A3)t/12]  + Ki~^— 

xJ  xy  U'X' 

— [(2A2A3  — A3)a;3i  + (A2  — 2A2A3)xi2]  + Ki~^- 


dp2  dp3 
dx  dx 
dp2  dps 
dy  dy 


It  now  follows  from  (3. 5, 3. 6)  that,  at  vertex  1, 


DXy  (2:1,  y\) 
DyX  Qy  (*^1  J 2/l) 


Dzy  + (Q 12  + ^12)2/31  + {Q31  ~ Gs\)yi2 

Dzx  + (Q12  + ^12)^13  + (Q3I  ~ C3i)x2l- 


In  view  of  Q12  + C12  = M21  — z2\,  Q3 1 — C31  = M3i  + z3 1 - by  the  definition  (2.7)  of  those 
quantities,  - a short  calculation  yields 


DXy  ( X\ , y\ ) 

DyX  {x\  1 y\) 


A^2i2/3i  + Mziyu  — Dxyzix, 


M2IX3I  + M3 \X\2  — DyXZly. 


This  - and  cyclic  substitution  (1.2)  - show  that  the  rHCT  function  z = z(x,  y)  indeed  assumes 
the  prescribed  derivative  values  at  the  vertices  of  the  footprint  triangle. 
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4.  The  Linear  Derivative  Condition 


The  salient  property  of  the  rHCT  element  - as  opposed  to  the  nonreduced  HCT  element  - is 
that  it  satisfies  (see  for  example  [10])  the  following 


(4.1)  Linear  Derivative  Condition : Along  each  edge  of  the  footprint  triangle,  the 
derivative  taken  in  the  direction  perpendicular  to  the  edge  varies  linearly  be- 
tween the  values  it  assumes  at  the  ends  of  the  edge. 


The  purpose  of  the  Linear  Derivative  Condition  is  to  ensure  that  adjacent  triangles  in 
a triangulation  fit  together  smoothly.  Indeed,  all  directional  derivatives  at  a vertex  can  be 
determined  from  the  - given  - gradient  at  that  vertex.  At  the  endpoints  of  the  common  edge, 
both  elevations  and  derivatives  in  the  direction  of  the  edge  are  therefore  predetermined.  Since 
the  rHCT  function  2 = z(x,y)  induces  a cubic  function  on  the  edge,  that  function  is  fully 
specified  - say,  by  Hermite’s  formula  - and  the  specification  is  the  same  in  both  the  adjacent 
triangles,  because  that  specification  is  solely  based  on  quantities  given  at  the  two  vertices 
of  the  boundary  edge.  Continuity  across  the  edge  therefore  holds,  as  does  commonality  of 
the  derivative  in  the  direction  of  the  edge.  In  order  to  ensure  smoothness,  that  is,  a common 
gradient  in  both  triangles  along  the  edge,  it  therefore  suffices  that  the  derivatives  perpendicular 
to  the  edge  are  also  fully  specified  by  their  values  at  the  vertices  of  the  edge. 

Now  if  the  perpendicular  derivative  along  an  edge  is  a linear  function,  that  is,  if  the 
Linear  Derivative  Condition  is  satisfied,  then  that  derivative  is  fully  specified  by  its  values 
at  the  vertices  of  the  edge.  This  is  all  that  is  needed  to  ensure  smoothness  across  triangle 
boundaries. 

A directional  derivative  of  2 = z(x,y)  perpendicular  to  the  edge,  say,  from  vertex  i to 
vertex  j is  given  by 


(4.2) 


dz  dz 

Viidi  ~ Xj% 


because  the  vector  (yjk,—Xjk)  is  perpendicular  to  the  direction  (xjk,yjk)  of  the  edge.  Any 
other  definition  of  the  perpendicular  derivative  differs  by  a constant  factor;  so  it  does  not 
matter  which  definition  we  choose. 

In  order  to  verify  that  the  rHCT  indeed  satisfies  the  Linear  Derivative  Condition  that  (4.2) 
is  a linear  function,  we  need  to  consider  only  the  three  functions  Vi  (2.8),  since  they  contain 
all  the  cubic  terms  of  the  function  2 = z(x,y).  In  particular,  we  need  to  consider  only  Vi, 
because  the  argument  extends  by  cyclic  substitution  (1.2)  to  the  two  remaining  functions  Vt. 

We  first  evaluate 


(4.3) 


( m dvl  \ 
xy  \V2Z^~X23^) 


3 


in  subtriangle  B\,  for  the  edge  connecting  vertices  2 and  3.  Here  X\  = 0,  whence  A2  + A3  = 1. 
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Since  the  gradients  of  p2  and  p3  vanish  along  that  edge  (see  3.4),  and  since  the  gradient  of  pi 
becomes 


AoA 


2/23  %23 


2/v3 


D ’ D 

-Lyxy  ^yx 


we  find  for  (4.3),  taking  into  account  Dyx  = —Dxy, 


(4.4) 

A3(2A2  - A3)(x23^3i  + 2/232/31 ) — X2(2X3  - A2)(o:232:i2  + 2/232/12)  + A2A3ATi(x23  + 1&)  ■ 
From  definitions  (2.7),  Ki(x\z  + 2/I3)  = 3(1/2  - T3).  Moreover,  the  scalar  products, 

S2  '■=  X23X12  + 2/232/12,  S3  :=  £31223  + 2/312/23? 
can  be  used  to  express  squares  of  edge  lengths  Lf: 

L\-L\  = 53-52. 

Expression  (4.4)  thus  reduces  to 


( A2  + A3)(A252  — X3S3) 


A2*S2  — X iS. 


3^3: 


and  is  therefore  linear  in  x and  y. 

Along  the  edge  in  subtriangle  B2,  we  have  A2  = 0,  A3  -f  Xi  = 1.  Only  the  gradient  of  p2 
does  not  vanish.  Similarly  to  (4.3)  we  find 

(A3  + Ai)(— X3LI)  = —X3L2, 

which  again  establishes  linearity.  In  subtriangle  jB3,  A3  = 0,  Ai  + A2  = 1.  Thus 

(Ax  + A2)(— X2LI)  = — A2T3. 

The  rHCT  is  therefore  seen  to  satisfy  the  Linear  Derivative  Condition. 


5.  Second  Derivatives  of  the  rHCT  Element 

We  now  aim  to  express  the  second  partial  derivatives 

d2z  d2z  d2z 
dx 2 ’ dx  dy  dy 2 
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as  linear  functions  of  the  barycentric  coordinates  Al5  A2,  A3.  We  start  with  the  correction 
functions  pl:  based  on  the  following  applications  of  the  Chain  Rule. 


d2pi  _ fdpi \ _d_/dpi\  dA2  _d_  fdpA  dXs 

dx 2 d\i  \ dx  J dx  + <9A2  \ dx  J dx  ^ dX3  y dx  J dx 

d2pi  _ d fdpi\  d ( dpi \ dX2  d f&Pi\  ^A3 

dxdy  dXi  \ dx  J dy  + dX2  \ dx  J dy  + <9A3  \ dx  J dy 

d2pi  d__  l dfh  \ _ 5Ai  _d_  f dpA  dX^  _d_  ( dpA  <9As 

dxy  dXi  V dy  ) - dy  dX2  \ dy  J dy  dX3  \dy  J dy 


From  (3.2,  3.3), 


(5.1) 


(5Ai  - 1 )y23+ 

d2  pi 

1 

2Aiy3i2/i2  + 2A  22/122/23  + 2A3?/23?/31 

in  B 

dx2 

' D2 

xy 

(A3  — A2)j/|1  + 2A22/3iyi2 

in  B 

, (A2  — A z)y\2  + 2A32/3i2/i2 

in  B 

(5.2) 


d2Pl 

dxdy 


1 

DxyD 


(5Ai  — 1)2/23^23  + 

X\(ysiXi2  + x^iVn)  + ^2(^12^23  + £122/23)  + A3(y23:r3i  + x23y3 1) 
(A3  — A2)y3i£3i  + X2(ysiXi2  + x^iyu) 


in  Bi 
in  B2 


K (^2  ~ ^3)2/12^12  + M(yzix\2  + £312/12) 


in  B3 


(5Ai  — 1)£23+ 


(5.3) 

d2pi 

dy2 

1 

" D2 

yx 

2Ai^:3i2:i2  + 2A2:ri2:r23  + 2X3x23x3\ 

(A3  — A2)  jTgj  + 2A2a;3ia;12 

in  B\ 

in  B2 

k ( A2  — A3)xJ2  + 2A3x3iX12 

in  B3 

Formulas  for  the  second  derivatives  of  the  two  remaining  correction  functions  p2  and  p3  follow 
by  cyclic  substitution  (1.2).  It  is  readily  verified  that  the  second  derivatives  are  continuous 
at  the  barycenter,  that  is,  each  subcubic  yields  the  same  value  for  Xi  = X2  = A3  = 1/3.  That 
continuity,  however,  does  not  generally  extend  along  every  common  edge  of  two  subtriangles. 
The  same  continuity  pattern  holds  for  the  full  element  (2.10). 
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We  now  turn  to  the  functions  V]  defined  in  (2.8)  whose  first  derivatives  have  been  calculated 
in  Section  3 (3.7).  It  follows  by  straightforward  applications  of  the  Chain  Rule  that 


(5.4) 


d2\\ 
dx 2 

d2Vx 
dx  dy 


d2VY 
dy 2 


D 2 

xy 


[A32/31  + (A2  — A3)t/12t/  31  — A2yi2]  + K\ 


d2  Pi  d2  p2  d2p3 
dx2  dx2  dx2 


Id  xyD  yx 


[d3yi3X3i  + (A2  — d3)(yi2X3l  + X122/31)  — A2?/21£  12  + 


K d<2pi  - d2p 2 + d2pz 
1 dxdy  dxdy  dxdy 


[A3X31  + (A2  — A3)xi2X3i  — A2rrf2]  + Hi 


D2 

yx 


d2pi  _ d2 p2  d2 p3 
dy2  dy2  dy2 


The  second  derivatives  of  the  two  remaining  functions  V2  and  V3  follow  from  the  above  by 
cyclic  substitution  (1.2).  We  can  now  express  the  second  derivatives  of  the  rHCT  function 
z = z(x,  y)  in  terms  of  the  second  derivatives  of  Vi:  i — 1,  2,  3. 


(5.5) 


d2z 

dx2 


d2z 

dxdy 


d2± 

dy: 


D2 

xy 


[0232/122/31  + ^31 2/232/12  + Ql22/3l2/23] 


+c„S  + c„S  + c,/ni 


dx2 


dx2 


dx2 


D xyD  yx 


[023(2/12^31  + £122/31)  + O31  (2/23^12  + £232/12)  + 012(2/31^23  + £312/23)] 


+c„|!li  + c„|^+c  *v> 


dxdy 


dxdy 


dxdy 


D2 

yx 


[023^12^31  + 031^23^12  + 0l2^31  ^23] 


d2V,  d2V2  d2\\ 

~ ~ t>3i  V - + Ui2' 


dy1 


dy: 


dy 5 


6.  The  Barycentric  Integrals 

The  functions  which  form  the  integrand  of  the  energy  integral  (1.1)  have  been  expressed  in 
terms  of  the  barycentric  coordinates,  but  the  variables  of  the  integration  are  still  x and  y. 
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More  precisely,  that  integrand  will  be  a quadratic  polynomial  in  Ai,  A2,  A3.  The  same  holds 
for  additive  components  of  the  integrand,  if  these  were  to  be  evaluated  separately. 

In  view  of  the  relation  Ai  + A2  + A3  = 1,  various  normalizations  of  such  polynomials  are 
possible.  For  instance,  the  polynomial  can  be  required  to  be  a homogeneous  quadratic  form, 
that  is,  contain  only  terms  of  degree  2.  Or  one  of  the  barycentric  coordinates,  say,  A3,  can  be 
expressed  in  terms  of  the  remaining  ones,  yielding  a nonhomogeneous  quadratic  polynomial 
in,  say,  X\  and  A2. 

Similarly,  the  cartesian  coordinates  x,  y can  be  expressed  linearly  in  any  two  of  the  barycen- 
tric coordinates.  For  instance,  eliminating  A3  in  relations  (2.1)  will  yield  such  an  expression  - 
and  its  cyclic  permutations  (1.2): 


x 

y 


x 

y 


X 

y 


( Xs  ) + ( Xl3  X23 

V 2/3  J { y 13  2/23 

f X1  \ ( X21  X31 

\ 2/l  / V 2/21  2/31 

( x2  \ + ( x32  X12 

v 2/2  ) v y 32 


Ai 

A2 

A2 

A3 

A3 

Ai 


Note  further,  that  the  determinants  of  the  above  linear  transformations  agree  with  the  previ- 
ously encountered  determinant  Dxy  (2.2)  of  the  linear  system  (2.1): 


det 

X13 

^23 

= det 

X21 

X31 

= det 

^32 

Xl2 

2/13 

2/23 

10 

2/31 

2/32 

2/12 

The  area  element  dx  dy  is  transformed  under  a linear  transformation 


into  the  area  element  det (T)dxdy.  We  thus  have 


(6.1) 


dx  dy  = DxydXidX2  = DxydX2dX3  — DxydX3dXi. 


In  view  of  A3  = 1 — Xi  — A2,  the  definition  (2.4)  of  the  subtriangle  B\  can  be  altered.  Up  to 
differences  of  measure  0,  we  have 


B\  { ( Ax , A2,  A3)  : 0 < Ai  < 1/3,  X\  < A2  < 1 — 2 Ax } , 


suggesting 

have 


limits  of  integration  over  those  subtriangles.  Indeed,  for  any  function  /(A1?  A2),  we 


/ /(Al5  A2)dA2  dXi  = /1_2Al  /(A1?  X2)dX2 

JB 1 JO  J Ai 


dXi 
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Using  this  procedure,  we  find  for  the  following  six  key  integrals: 


'B  i 


1 d\2 d\\  — 


'Bl 


\\dX2  d\\  — 


\2d\2  d\\  — 


'Bx 


f \\d\2  d\i  = 

B 1 

\1\2dX2  d\\  — 


,B1 


A9C/A2  dX\  — 


r l-2Ai 
’Xi 

rl-2Xi 

Ai 

rl-2X1 

'A: 

•1— 2Ai 


1 dX2 

X\dX2 

A2C/A2 


dXi  = t 
6 


dXl  = 


54 


dX1  = — 
27 


[ Xfd\2 

J Ai 

r\—2Xi 


dXl  = 


324 


Aj 
rl—2Xi 

X-i 


XiX2dX2 


dX1  = 


648 


X\dX2 


dX-i  — 


13 

324 


These  key  integrals  and  their  cyclic  equivalents  are  all  that  is  needed  for  the  energy  formula 
to  be  derived  in  the  next  section. 

In  general  and  for  checking  purposes,  it  may  be  useful  to  provide  a complete  list  of  barycen- 
tric  integrals  in  subtriangles  Bt.  All  the  remaining  integrals  can  be  derived  from  the  above 
key  integrals  using  the  relations  Ai  + A2  + A3  and 

(6.2)  dXi  + dX2  T dA3  — 0. 

Since  A3  = 1 — Ax  — A2, 

/ A3dA2  dXi  — / \dX2  dX\  — / X\dX2  dX\  — / X2dX2  dX\ 

J B\  J B\  J B\  J B\ 

11  2 _ 2 
6 _ 54  _ 27  _ 27 

Analogous  calculations  yield  the  remaining  integrals  over  subtriangle  B\  based  on  the  area 
element  dX2dXi.  Integrals  over  the  subtriangles  B2  and  B3  follow  from  integrals  over  B\  by 
cyclic  substitution  (1.2). 

Suppose  indices  z,  j,  k are  in  cyclical  order.  Then  substituting  — dX3  — dX^  for  dXi  yields 


dXj  dXi  — — dXj  dXj  — dXj  dX^. 

Since  area  elements  change  orientation  when  the  sequence  of  the  differentials  is  switched,  and 
since  consequently  also  dX3  dX3  = 0,  it  follows  that 

dX2  dX\  — d\3  dX2  — dX\  dX$. 
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As  a result,  it  is  immaterial  which  of  the  three  area  elements  is  specified,  so  that  only  the 
integrand  is  of  interest.  This  permits  the  display  of  our  barycentric  integrals  in  Table  (6.3). 

(6.3)  Barycentric  integrals: 


B1 

b2 

S3 

B 

1 

1 

1 

1 

1 

6 

6 

6 

2 

Al 

1 

2 

2 

1 

54 

27 

27 

6 

^2 

2 

1 

2 

1 

27 

54 

27 

6 

^3 

2 

2 

1 

1 

27 

27 

54 

6 

1 

13 

13 

1 

324 

324 

324 

12 

13 

1 

13 

1 

324 

324 

324 

12 

A! 

13 

13 

1 

1 

324 

324 

324 

12 

V 

v 

to 

5 

'648 

5 

648 

17 

648 

1 

24 

A1A3 

5 

648 

17 

648 

5 

648 

1 

24 

A2A3 

17 

648 

5 

648 

5 

648 

1 

24 
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Additional  symmetries  are  apparent.  In  subtriangle  Bx , the  roles  of  variables  X3  and  Xk  are 
interchangeable  as  far  as  integration  is  concerned.  That  is  at  the  root  of  relations 


f f{Xl,Xj)dXjdXi 

J Bx 


j f (A^  Xk) dXj  dXi, 
J Bi 


j,k  ^ i. 


B 


f(Xu  X2)dX2  dXi  = 


•1  — Ai 


f{X1,X2)dX2 


dAi 


Those  results  are  also  displayed  in  Table  (6.3).  They  should  be  the  sum  of  the  integrals  of  the 
same  integrand  over  the  three  subtriangles  Bx. 


7.  Derivation  of  the  rHCT  Energy  Formula 


In  this  section,  we  will  carry  out  the  integration  prescribed  by  the  formula  (1.1)  for  the 
surface  energy  of  an  almost  flat  thin  plate  for  a single  triangular  rHTC  element.  That  element 
is  determined,  as  set  forth  in  the  previous  sections,  by  the  vertex  data  supplied  at  the  corners 
of  the  given  triangle  B. 

Let  B = B0U BiU B2U B3  be  the  given  triangle  with  its  barycentric  partition.  The  surface 
energy  E of  the  surface  element  is  then  given  by 


(7.1) 


E := 


dx  dy  , 


where  z = z(x.y)  is  the  rHCT  function  (2.10).  Our  goal  is  to  express  E in  closed  form  in 
terms  of  the  elevations  and  gradients  at  the  triangle  vertices  as  well  as  the  triangle  geometry. 
Energy  E is  the  sum  of  the  energies  Ex  in  the  three  subtriangles  Blx  i — 1,  2,  3, 


E — Ei  + E2  + £3, 


and  will  be  determined  in  that  fashion.  Once  an  expression  for  E\  has  been  obtained,  energy 
expressions  for  E2  and  E3  follow  by  cyclic  substitution. 

The  second  derivative  expressions  (5.1)— (5.5)  for  the  pi  s,  V^s,  and  hence  z,  are  linear  in 
Ai,  X2  and  A3.  Also,  since  Xi  + X2  X3  = 1,  we  can  replace  A3  everywhere  with  1 — Ai  — A2, 
resulting  in  second  derivative  expressions  as  follows  which  are  linear  in  Ai  and  A2. 
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1 


d2z 
dx 2 


D2 

xy 


(Cqxx  “1“  “I-  C2XX-A2) 


d2z 
dx  dy 


T-,9  (O) xy  T ^lxyAj  + C2xyX2) 

Uxy 


JJ2  (C°yy  C!yy/^1  C2yy^2)  ■ 

Uxy 

The  coefficients  c0xx , cixx, ...  of  the  above  linear  form  refer  to  subtriangle  Bi:  and  are,  of  course, 
different  for  the  other  subtriangles.  Bi  of  the  barycentric  partition. 

The  second  derivative  expressions  also  show  that  the  coefficients  c...  themselves  are  homo- 
geneous quadratic  forms  in  the  coordinate  differences  all  of  which  can  be  expressed 

linearly  in  x2i,  X13,  yn,  2/31  since  :r32  = — rr13  — x2i,  etc.  Due  to  the  structure  of  the  deriva- 
tive expressions  only  nine  coefficients  ai3k  occur  as  coefficients  of  those  quadratic  forms.  For 
i = 1,2,3, 


d2z 
dy 2 


where 


Qxx  — ai22Vi2  + 2 dmynVzl  + ai33Vl\ 

CiXy  = CLi22X2iy\2  + Ui23(^2l2/31  + #132/12)  + a033#13j/31 

2 2 
ciyy  — &i22#21  + 2<2j23#21#13  + az33#i3  , 


&022  = 

= —2  Qzi  + 2Ci2 

— 

(K 1 + 1)C23  -f  (K2  — 7)C3i 

&023  = 

= — Ql2  + ^23  ~ 

Q31  + 3(7* 

12  - (2A 

1 + 3)Q 

>3  + (K2  ~ 

a033  r 

- ~ 2Qi2  + 4(7*12 

— 

(3/fi  - 

l)(7*23  + 

(AT2  - 5)C31 

&122  - 

= — {Ks  + 9)Ci2 

+ 

(bKi  + 

3)  C23  — 

(4A-2- 

18)(7*3i 

&123  ~ 

= — (2iT3  + 12)C 

12 

+ (7  K, 

+ 3)C*23 

- (3A'2 

— 15)(7*3i 

C-133  = 

- — (3j F73  + 15)<7 

12 

+ (7  Kx 

3)  (7*23 

- (2A-2 

- 12)C*31 

&222  = 

= + (-^3  + 3)Ci2 

— 

2KyC2Z 

(#2  - 

■ 3)C*3i 

^223  = 

r + 3)Ci2 

+ 

6O23  — 

(AT2  - 3 )C31 

&233  “ 

= +(A3  + 3)Ci2 

+ 

2K\C23 

-(^2- 

■ 3)  C31 

Note  that  the  coefficients  are  homogeneous  linear  forms  in  the  quantities  whose  co- 

efficients depend  only  on  geometric  information,  K\  — 3(L\  — L\)/L\, ...,  independent  of  the 
choice  of  the  coordinate  system. 

The  coefficients  a...  above,  as  well  as  the  subsequent  results  reported  here,  were  mostly 
determined  with  the  help  of  Mathematica  [16]  and  verified  by  comparison  with  hand  calcula- 
tions. 
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So,  the  integrand  in  the  energy  integral  (7.1)  can  be  expressed  over  subtriangle  B\  in 
barycentric  coordinates  as  a quadratic  polynomial  in  the  barycentric  coordinates  Ai,A2: 

(ePzV  ( &z  \2  fd^zV 

\dx2 ) + \dxdyj  \9y2/ 

= J^-{CQXX  + C1  xx Ai  “1“  C2xX^2^  4“ 

^Xy 

2 

p.^  (O) Xy  4~  C\Xy\\  4“  C2XJ/A2)  4~ 

UXy 

(C0 yy  4~  C\yy\\  4*  C2yyX2) 

Uxy 

+ 2g0iAi  4-  2^02^2  + Qn^i  4-  2^12A1A2  4-  922A2)  • 

^Xy 

The  coefficients 


Qij 


CixxCjxx  d-  2CixyCjxy  4-  C-iyyCjyy  ■> 


0 < i < j < 2, 


of  that  quadratic  function  are  themselves  quadratic  forms  in  the  nine  coefficients  a...  with 
coefficients  that  are  homogeneous  biquadratic  forms  in  coordinate  differences.  To  prepare  for 
that  evaluation,  we  calculate  the  three  products  in  the  above  formula  separately. 

CiXXCJXX  = 0'i22^j22Vi2  4-  ^i22aj232yl2y3l  4"  <*£22 aj33 2/12^31  + 

^i23aj222y\2y3l  + ^i23aj23^yi2yll  + ai23^j332yi2yh  + 
a i33aj22yi2yll  + ai33aj232ynyl\  + 0*33  Oj'33  2/31 


C-ixyCjxy  ^i22^j22x\iV\2  + ai22&]23  (^ll  2/122/31  + ^21  ^ 13  y 12  ) az22  <^y33  ^2 1 ^13  2/12  2/31  + 

0'i23®jj22  (^21 2/122/31  + ^21^132/12)  ”b 

^23^23(^212/31  4-  2^21^132/122/31  + ^132/12)  4" 

&i23&j33 (^21^132/31  4*  ^132/122/31 ) 

ai33aj22x21x13yi2y3l  4"  Gf33%23  (^21^132/31  + ^132/122/31 ) + ai33^j33x\zy\i 


ciyycjyy  = ai22aj22x2l  4"  ai22aj23'2x2ixl3  + ai22aj33x21X13  4~ 

^£23  ^'22  2^21*^13  + ^£23^2343:21^13  + G£23^j3323:2l3;j3  + 
tti33aj22x21x13  4*  &i33aj23 2^21^13  4*  ^£33^33^13  • 
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By  substituting  into  the  additive  expression  for  coefficient  and  collecting  terms, 

Qij  = Gi22«j22  (^21  + Vn)2  + 

<L22%23  2(^21  + 2/12) (X2\X\3  + 2/l22/3l)  + 

^22^33(^21^13  + 2/122/31  )2  + 

&z23aj22  2(^2!  + 2/12)  (^21^13  + 2/122/31 ) + 


Since 


we  find 


&i23aj23  2 (£21X13  + 2/122/31  )2  + (^21  + 2/12)  (X13  + 2/fl ) 
&i23aj33  2(Xj3  + 2/31)  (^21^13  + 2/122/31 ) + 

&i33®'j22  (x21x13  + 2/122/31  )2  + 

tii33aj23  2(^13  + 2/31 ) (X2l£l3  + 2/122/31 ) + 

flz33%33(^i3  + 2/fl) 2 " 


+ 


Qij 


£2l£l3  + 2/122/31  — 


L?  - L|  - Ll 


— &i22aj22  (£3)  + 

2ai22dj23 


r 2 r 2 t 2 

^1  ^2  L3  1 r 2 


Q'i22&j33 
2az23aj22 
2ai23Clj23 
2dl23Clj33 
&i33aj22  I 
2cii33aj23 


L\-L\  - L\ 


Li  + 


+ 


L\-Ll-L\,  t2 


Li  + 


T 2 _ r2  _ r 2 \ 2 
^1  -^2  ^3  \ 


+ L9L 


2 r 2 


2^3 


'£?  - ^ - £§' 


7-2  _ r2_  r2^ 
^1  ^2  -^3 

2 

fL\-Ll-L\ 


L\  + 


+ 


Z|  + 


ai33&j33  (jA ) • 


+ 
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Note  that  the  coefficients  q„  thus  represent  purely  geometric  information  independent  of  the 
choice  of  coordinate  system,  because  the  coefficients  a...,  too,  are  independent  of  the  choice 
of  the  coordinate  system,  as  seen  earlier.  This  indicates  that  the  expression  for  the  surface 
energy  is  independent  of  the  choice  of  the  coordinate  system,  as  we  would  expect.  (Note  also 
that  the  g..’s  will  be  different  for  each  subtriangle  Bx  in  the  partition  of  triangle  £,  since  the 
coefficients  c ..  are.) 

It  remains  to  substitute  the  values  of  the  key  barycentric  integrals.  So,  the  energy  integral 
over  subtriangle  B\  is 

1/1  1 4 1 5 13  \ 

El  ~ Dlv  U9o°  + 27 901  + 27 902  + 3249u  + 32i912  + 3 2 4 922J  ’ 

xy 

where  the  qj s are  those  appropriate  for  B\.  This,  finally,  results  in  a homogeneous  quadratic 
form  in  the  derivative  quantities  QlJ:  CXJ 

E\  = ( g\Cl  ex  023023  T #lcx  02  023031  T #1ciC3  023012 

T 9lc\q\  023023  + 9lciq2  023031  + #10193  023012 

+ #lc2c2  031031  + #lc2c3  031012  + #lc29x  031023 

T 9lc2q2  031031  T #10293  031012  T #10303  012012 
+ #lc39i  0L2023  + #10392  012031  + ^10393^12012 
+ ^19x9x023023  + #19x92023031  + #19x93  023012 
+ #19292  031031  + #19293031012  + #19393012012  )/ D\y 

with  coefficients  that  are  homogeneous  rational  functions  - as  displayed  below  - in  the 
squares  L2  of  the  side  lengths  of  triangle  B.  The  coefficients  of  these  rational  functions  are 
rational  numbers. 

#ioici  = (3L^-6L^L2  + 88Lt^-12L2^  + llL^-6L5^-140Tt^2^3  + 12^L2L3 

- 40L62L2  + 88L1L3  + 12 L\L\L\  + 58 L\L\  - \2L\L%  - 40 L\L\  + 11L|)/(36LJ) 

#ioio2  = {nL\  + ^L\Ll  + ZmL\L\-^L\Ll  + hL\-^L\L\-Z^L\LlL\ 

+ 72 L\L\L\  - 28 L\L\  + lOOLjLj  + 38 L\L\  - 24 L\h\  - 12 L\L\  - ZL\)  / (72L\L\) 

9io  103  = (11L?  - 84 L\L\  + 100 L\L\  - 24 L\L\  - 3 L\  + 68LflJ  - 3 ^L\L\h\  - 12 L\L\ 

+ 300 L\L\  + 72L\l?2L\  + 38L42L$  - 4 8L\h\  - 2 8L\h\  + $L\)  / {72L\L\) 
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9\clql 


9lclq2 


9lclq3 


9lc2c2 


9lc2c3 


9lc2ql 


9lc2q2 


9lc2q3 


9lc3c3 


9lc3ql 


9lc3q2 


(-7 L\L\  - 6 L\L\  + 13 L\  + 7 L\L\  - 31  L\L\  + 6 L\L\ 

+ 31  L\L\  - 13Z|)/(18£?) 

(-2 L\  + 23 L\L\  - 20 L\L\  — L%  — 15 L\L\  + 40 L\L22Ll 

- L\L\  - 20 L\L\  + 5 L\L\  - 8L\)/(18L\) 

(2 L\  + 15 L\L\  + 20 L\L\  + 3 L%  - 23 L\L\  - AQL\LlLl 

- 5 L\L\  + 20 L\L\  + L\L\  + Z,|)/(18Z,|) 

(131®  + 78 L\L\  + 124 L\L\  - 50 L\L\  + 3 L\  - 26 L\L\  - 92 L\L\L\ 

+ 62 L\L\L\  - 4 L\L\  + 14 L\L\  + UL\L\L\  - 2 L\L\  + l|)/(72L|) 

(-6L?  - 18 L\L\  + 39 L\L\  - 16 L\L\  + L\  - 18 L\L\  - 168 L\L\L\  + 28 L\L%L\ 
+ 2 L\L\  + 89L\L\  + 28L\L2L\  - 6 L42L43  - 16 L\L%  + 2 L\L\  + L?)/(36Z|l|) 

(-2 L\  - 11  L\L\  + 8 L\L\  + 5 L\  + 4 L\L\  - 6 L\L%L\ 

- 22 L\L\  - 2 L\L\  + 17l|l4)/(18Z|) 

(5 L\  + 12L\L\  - 19 L\L\  + 2 L%  - 5 L\L\  + 32 L2L22L23 

- 3 L\L\  - L\L\  + Ll)/{\8Ll) 

(3 L\  + 28 L\L\  + 11  L\L\  - 2 L\  - 7 L\L\  - 28L\L\L\ 

+3 L\L\  + 5 L\L\  - Lt)l(\8Ll) 

(131?  - 26 L\L22  + 14 L\L\  - 2 L\L%  + L%  + 18L\l\  - 92 L\L\L\  + UL\L\L\ 

+ 124 L\L\  + 62 L\L%L\  - 50 L\L\  - 4 L\L%  + 8L\)j{12L\) 

(2 L\  - 4 L\L\  + 2 L\L\  + 11  L\L\  + 6 L\LlL\  - 17L\L23 
-8L\L\  + 22L2L\-bLl)/{l8Ll) 

(-3 L\  + 7 L\L\  - 5 L\L\  + L%-  28 L\L\  + 28 L\L22L\ 

- UL\Ll-8L2L\  + 2Ll)/{l8Ll) 
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<?1C3,3  = (-5  L\  + 5 L\L\  + L\L\  -L\-  12  L\L\  - 32  L\L\L\ 

+ 1§L\L\  + ZL\L\  - 2Lf)/(18Lf) 

= (L\-2L\Ll  + L\-2L\Ll  + 6LlLl  + Li)/l2 

91.1.2  = {-L\  + 2L\Ll-L\-2L\Ll-2L%Ll  + ZL\)ie, 

91.1.3  = {-L\-2L\Ll  + ZL\  + 2L\Ll-2LlLl-L\)ld 

91.2.2  = {L\-2L\Ll  + L\  + e,L\Ll-2LlLl  + L\)/l2 

91.2.3  = {ZL\-2L\L\-L\-2L\L\  + 2L\L\-L\)IZ 

91.3.3  = {L\  + SL\Ll  + L\-2L\Ll-2LlLl  + L\)l\2 

The  above  coefficients  reflect  inherent  symmetries  with  respect  to  interchanging  the 
indices  2 and  3.  Some  coefficients  are  invariant  with  respect  to  this  transposition.  Others 
become  the  positive  or  negative  of  the  coefficient  whose  subscripts  have  been  interchanged. 
The  negative  outcome  occurs  for  the  coefficients  associated  with  mixed  products  CijQkh  and 
is  due  to  the  fact  that  the  derivative  quantities  CZJ  switch  signs  as  indices  i and  j are  inter- 
changed, whereas  the  derivative  quantities  Qki  remain  unchanged  under  index  transposition. 
The  following  transitions  are  generated  by  the  interchange  2 *->•  3: 

(7.2)  <?lclcl  +9lclch  9lclc2  — > +giclc3>  9lc2c2  +#lc3c3:  <?1c2c3  “ > +<?1c2c3: 

9lclql  —9lclqli  9\c\q2  ~ 9lclqS,  9lc2ql  ~9\c3ql j 9\c2q2  ~ > ~ 9\c3q3i 

9lqlql  ~^~9lqlql  i 9lqlq2  +9lqlq3:  9lq2q2  +9lq3q3i  9lq2q3  +9lq2q3  ■ 

The  remaining  energy  expressions  E2  and  E$  follow  from  the  above  one  for  E\  by  cyclic 
substitution.  Adding  those  three  yields  the  complete  energy  expression 

E — ( 9C1C1C23C23  + 9cic2 C23C31  + 9C1C3C23C12 

+ 9c\q\  C23Q23  + 9ciq2E'2^Qz\  + 9cxq3  C23Q12 
+ 9c2c2  C31C31  + 9c2C3C3lCl2  + Pc2<7i(-3lQ23 

+ ^0292^31^31  + 9c2qz^3lQl2  + 9c3c3C  12C12 

+ 9c3qiCi2Q23  + 9czq2CuQ3\  + 9c3q3CnQ  12 
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+ 9qiqiQ2:iQn  + 9qiqqQxiQ‘‘A  + SqiqqQ'l'lQl’l 
+ S9292Q31Q3I  + 9q2q3Q3lQl2  + 9q3q3Ql2Q\2  ) / &ly  , 

where 

gcici  = (2L\-llL\Ll  + WL\L\  + $L\Ll  + SLl-UL\Ll-2§L\LlLl-$L\L\Ll 

- 1 8L62Lj  + 50 L\L\  - 9L\L22Ll  + 24 L\L\  + 9 L\L\  - 18 L\L\  + 6Lf)/(12Z,?) 

gClC2  = (3 L{  + 4 L\L\  + 98L\L\  + 4 L\L\  + 3 L\  - 24L\L\  - 40 L\LlLl  - 40 L\L\L\ 

- 24 L%L\  + 36 L\L\  - 56 L\L\L\  + 36 L42L43  - 12 L\L\  - 12 L\L%  - 3L|), / (12L\L\) 

gClC3  = (3 L\  - 24 L\L22  + 36Z?L|  - 12 L\L\  - 3 L\  + 4L\L\  - 40 L\L22Ll  - 56 L\L\L\ 

- 12 L\L\  + 98 L\L\  - 40 L\L%L\  + 36 L\L\  + 4 L\L\  - 24 L\L\  + ZL\) / {12L\L\) 

9cun  = {L\L\-2L\L\  + L\-L\L\-ZL\L\  + 2L\L\  + %L\L\-L%)I(2L\) 

3c„2  = {-L\  + 2L\Ll-Lt  + 2L\Ll  + 2LlLl-L\)/2 
gciq3  = (L\-2L\Ll  + L\-2L\Ll-2LlLl  + L\)/2 

gcm  = {8L\  + 9L\L\  + h9L\L\~\lL\L\  + 2L\-\8L\L\-9L\L\L\-28L\L\L\ 

- 11  L\L\  + 24 L\L\  - 9 L\L\L\  + 50 La2L\  - 18 L\L\  + 9 L\L%  + 6L|)/(12L^) 

5C2C3  = (-31?  - 12 L\L\  + 36 L\L\  - 24 L\L%  + 3 L\  - 12L?Lj?  - 56 L\L\L\  - 40 L\L\L\ 

4-  4 L\L\  + 36 L\L\  - 40 L\L\L\  + 98i^L£  - 2 4L\Ll  + 4 L\L%  + 3Zj)/(12l|L£) 

gcw  = (L4  - 2L\L\  + L\-  2L\L\  - 2L\L\  + L\)/2 


3c2,2  = {-L\  + 2L\L\-L\L\  + ZL\L\  + L\L\-ZL\L\-2L\L\  + L%)I{2L\) 
ffc293  = {-L\  + 2L\Ll-Ll  + 2L\Ll  + 2LlLl-Lt)/2 

gC3C3  = (61?  - 18 L\L\  + 24L\L\  - 18 L\L\  + 6 L\  + 9 L\L\  - 9 L\L22L\  - 9 L\L\L\ 

+ 9 L\L\  + 50 L\L\  - 26 L\L22L\  + 50i|l£  - 11  L\L\  - 11  L\L%  + 2I|)/(12L|) 
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9czqi 


9 C3<72 


9c$qz 


9q\q\ 


9qiq2 


9q\qz 


9qzqz 


9 <72  <73 


= (6 L\  - 1 8L?Z|  + 24 L\L\  - 18 L\L\  + 6l|  + 9 L\L\  - 9 L\L\L\  - 9 L\L\L\ 

+ 9 L\L\  + 50 L\L\  - 26 L\L%L\  + 50 L\L\  - 11  L\L\  - 11  L%L%  + 2L\)/{12L\) 

= (L\  - 2 L\L\  + L\-  2L\L\  - 2 L\L\  + L\)/ 2 

= (L\  - 3 L\L\  + 3 L\L\  -I*-  2 L\L\  + 2L\L\  + L\L\  - L\L\)I(2L\) 

= (L?  - 2 L\L\  + L4  - 2 L\L\  + + L\)j 4 

= (-L4  + 2i?I^  - L\-  2L\L\  - 2L\L\  + 3L|)/2 

= (-Li  - 2 L\Ll  + 3 L\  + 2L\L\  - 2 L\L\  - £4)/2 

= (if  - 2L?L|  + -^2  + 61?^  - 2 L\L\  + L43)/4 

= (3  Li  - 2L\L\  -Li-  2L\Ll  + 2 L\L\  - L\)j2 


9, 3,3  = + f>L\L\  + L\  - 2L\L\  - 2L\L\  + 1*) /A 


Again  the  symmetries  of  the  problem  are  reflected  in  transformations  of  the  above  coeffi- 
cients under  permutations  of  the  indices  i = 1,2,3.  The  relationships  (7.2)  obtain  here,  too, 
for  the  index  interchange  2 3.  In  addition,  the  corresponding  relationships  are  observed 
for  the  other  two  transpositions  3 f>  1 and  1 <-»  2.  Those  transpositions  can  be  combined 
to  generate  all  other  permutations.  Note  that  all  coefficients  of  the  form  gcjqk,  j 7^  are 
symmetric  functions  in  the  squares  and  are  equal  up  to  their  sign. 
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